import math

#################################################
#Constantes                                     #
#################################################
# - pi
pi = math.pi
# - Permeabilidad magnetica en el vacio
mu0=4.0*pi*1.e-7  #Newtons/Amperio2
# - Constante de Boltzmann
kB = 1.3806504e-23  #J/K

#################################################
#Propiedades del campo externo                  #
#################################################
# - Campo magnetico
H0=10000  #Amperio/metro 

#################################################
#Propiedades de la disolucion                   #
#################################################
# - Viscosidad de la disolucion
vis_disolucion= 200.0e-3 #Pa*s
# - Permeabilidad relativa de la fase continua
mucr = 6.0 #adimensional
# - Temperatura
Temp = 273 + 25
# - kBT
kBT=kB*Temp
#kBT=1.0

#################################################
#Propiedades de las particulas                  #
#################################################
# - Tamano de particula
radio=400.0e-9          #metros
#radio=0.5
#radio=1.0
diametro=2.0*radio      #metros
# - Volumen de particula
vol_part= 4./3.*pi*math.pow(radio,3.0)  #m3
# - Densidad de las particulas
rho_part = 7.8e3 #kg/m^3
#rho_part = 1.0e3 #kg/m^3
# - Masa de la particula
m_part = vol_part*rho_part #kg
#m_part=1.0
# - Permitividad magnetica relativa de la particula
mupr=100.0 #adimensional
# - Factor de contraste
beta = (mupr-mucr)/(mupr+2.0*mucr) #adimensional
# - Momento magnetico de las particulas
m=4.0*math.pi*mu0*mucr*beta*math.pow(radio,3.0)*H0 #Newton*m2/Amperio (revisar)
print
print '**********PROPIEDADES DE PARTICULA***********'
print "Particle diameter (m):",diametro
print "Particle mass (kg):",m_part
print "Magnetic moment (N*m2)/A: ",m

#################################################
#Interaccion magnetica                          #
#################################################
# - Coeficiente de energia de interaccion U0
u0 = mu0*mucr*beta*beta*math.pow(radio,3.0)*H0*H0 #Julios
F0=3.*u0/radio #Newtons
print 
print '**********PROPIEDADES DEL CAMPO***********'
print 'u0 (J)',u0
print 'F0 (N):',F0

#################################################
#Ecuacion de Langevin                           #
#################################################
# m*a = -gamma*v +random + Forces
# D = kBT/gamma (Coeficiente de difusion)
# gamma = 6*Pi*etaMedio*radio(Coeficiente de friccion)
# - Coeficiente de friccion gamma
gamma = 6.0*pi*vis_disolucion*radio/m_part #1/s
# - Coeficiente de difusion
Dif= kBT/gamma #m2/s
print
print '**********PROPIEDADES DE ECUACION DE LANGEVIN***********'
print 'Coeficiente de friccion gamma (N*s/m):',gamma
print 'Coeficiente de difusion D (m2/s)',Dif


#################################################
#Unidades adimensionales                        #
#################################################

# - Unidad de longitud radio
Lr=radio
# - Unidad de masa
Lm=m_part
# - Unidad de tiempo
Lt =math.sqrt(m_part*radio/F0)
# - Gamma reducido
gamma_red=gamma/Lt
# - Difusion reducida
Diff_red=Dif*Lt/kBT
print "Gamma program: ",kBT/m_part/Dif*Lt
print "Difusion program",1./(kBT/m_part/Dif*Lt)
print "treal = tred*",Lt," segundos"
print "Coef. Friccion reducido: ",gamma_red
print "Coef. difusion reducido: ",Diff_red
print "u0 reducido",u0/kBT

gamma_red2=gamma*math.sqrt(radio/(F0*m_part))
print "gamma_red2: ",gamma_red2
print "Diff_red2: ",1./gamma_red2
print math.sqrt(m_part*radio/F0)*gamma

